function  [xhat,gamma,xmu,ivar,sqrtvar,var,eps,Y]=batchnorm_forward(x,gamma,beta,eps)
[N, D]=size(x);
% 1. calculate mean
mu=1./N*sum(x,1);
% 2. substract mean vector
xmu=bsxfun(@minus,x,mu);
% 3. denominator
sq=xmu.^2
% 4. calculate variance
var=1./N*sum(sq, 1);
% 5. add eps for numerical stability
sqrtvar=sqrt(var+eps);
% 6. invert sqrtvar
ivar=1./sqrtvar;
% 7. execute normalization
xhat=bsxfun(@times,xmu,ivar);
% 8. transformation steps
gammax=bsxfun(@times, gamma, xhat);
% 9.
Y=bsxfun(@plus, gammax, beta);
% 10. store intermediate
end

function [dx,dgamma,dbeta]=batchnorm_backward(dY,xhat,gamma,xmu,ivar,sqrtvar,var,eps)
% unfold
[N,D]=size(dY);
% step9
dbeta=sum(dY,1);
dgammax=dY;
% step8
dgamma=sum(bsxfun(@times,dgammax,xhat),1);
dxhat=bsxfun(@times,dgammax,gamma);
% step7
divar=sum(bsxfun(@times,dxhat,xmu),1);
dxmu1=bsxfun(@times,dxhat,ivar);
% step6
dsqrtvar=-1./(sqrtvar.^2).*divar;
% step5
dvar=0.5*1./sqrt(var+eps).*dsqrtvar;
% step4
dsq=1/N*bsxfun(@times,ones(N,D),dvar);
% step3
dxmu2=2*xmu.*dsq;
% step2
dx1=(dxmu1+dxmu2);
dmu=-1*sum(dxmu1+dxmu2);
% step1
dx2=1/N*bsxfun(@times,ones(N,D),dmu);
%  step 0
dx=dx1+dx2;
end